Astronomy & Astrophysics manuscript no. LeBorgne09 v2 


© ESO 2009 


July 2, 2009 





o 
o 

(N 



(N 



O 



> 
m 

00 

cn 

o 

O 



X 



Cosmic star-formation history 
from a non-parametric inversion of infrared galaxy counts 

D. Le Borgne^'^, D. Elbaz^ P. Ocvirk' ^ C. Pichon^ 
CEA/Saclay, DSM/IRFU/SAp, F-91191 Gif-sur-Yvette, France 

Institut d'Astrophysique de Paris, UMR7095 CNRS, UPMC, 98bis boulevard Arago, F-75014 Paris, France 
Astropliysikalisclies Institut Potsdam, An der Stemwarte 16, D-14482 Potsdam, Germany 



Received...; accepted... 

ABSTRACT 

Aims. This paper aims at providing new conservative constraints on the cosmic star-formation (SF) history from the empirical mod- 
eling of recent observations in the mid and far infrared. 

Methods. We present a new empirical method based on a non-parametric inversion technique. It primarily uses multi- wavelength 
galaxy counts in the infrared and sub-mm (15, 24, 70, 160, 850 fim), and it does not require any redshift information. This inversion 
can be considered as a "blind" search for all possible evolutions and shapes of the infrared luminosity function of galaxies, from 
which the evolution of the star-formation rate density (SFRD) and its uncertainties are derived. The cosmic infrared background 
(CIRB) measurements are used a posteriori to tighten the range of solutions. The inversion relies only on two hypotheses: (1) the lu- 
minosity function remains smooth both in redshift and luminosity, (2) a set of infrared spectral energy distributions (SEDs) of galaxies 
must be assumed, with a dependency on the total luminosity alone. 

Results. The range of SF histories recovered at low redshift is well-constrained and consistent with direct measurements from various 
redshift surveys. Redshift distributions are recovered without any input into the redshifts of the sources making the counts. A peak of 
the SFRD at z 2 is preferred, although higher redshifts are not excluded. We also demonstrate that galaxy counts at 160 fim present 
an excess around 20 mJy that is not consistent with counts at other wavelengths under the hypotheses cited above. Finally, we find 
good consistency between the observed evolution of the stellar mass density and the prediction from our model of SF history. 
Conclusions. Multi-wavelength counts and CIRB (both projected observations) alone, interpreted with a luminosity-dependent library 
of SEDs, contain enough information to recover the cosmic evolution of the infrared luminosity function of galaxies, and therefore 
the evolution of the SFRD, with quantifiable errors. Moreover, the inability of the inversion to model perfectly and simultaneously 
the multi-wavelength infrared counts implies either (i) the existence of a sub-population of colder galaxies, (ii) a larger dispersion of 
dust temperatures among local galaxies than expected, (iii) a redshift evolution of the infrared SED of galaxies. 

Key words. Galaxies: high-redshift- Galaxies: evolution - Galaxies; formation - Infrared: galaxies - Submillimeter - Galaxies: 
luminosity function 

1. Introduction used extensively to measure the star-formation rate (SFR) of 

both nearby and distant galaxies. The evolution of the infrared 
Some key questions remain concerning the formation of galax- (j^^ luminosity functions (hereafter LF), parameterized in shape 
les, such as when and how galaxies formed their stars over the g ^-^^^ Schechter functions), has been measured up to z = 2 



^ . past 1 3 Gyr. Thanks to recent ultra-deep surveys at various wave- ^-^^^^ ^^^-^^ ^ parameterization for the evolution that can be both 
- - . lengths, some phenomena are now quite accurately measured in lu minosity ((1 or in density ((1 From these stud- 

and described, at least at relatively low redshift. For instance, it is jgs dLe Floc'h et alJ l200l iBabbedge et alJ & ICaputi et al.l 



well-established that massive galaxies have experienced most of 12007D, several valines of and ao have been measured in var- 
their SF activity at early epochs, whereas the SF activity in small -^^^ ^^^^j^jfj ^^^^^^ ^^^^ ^^^-^^ jj^^ evolution of the 
galaxies keeps a more constant level, on average. This so-called gp^^ g^^j^ ^^^y.^ ^ ^^^-^ ^^^-^ understanding 



downsizing" has been subject to many studies over the past few g^j^^y evolution, but they are generally limited to relatively 
ears (Madau et al. 1996; Lilly et al. 1996; Steidel et al. 1999; ^g^^hifts for two reasons. First, they are very expensive in 



. ._ ^ \^\JlJJ.±iX LiJ L vy \_/ X \^ LI. I \ / I I.I , i iX iJ L, l.ll\^Jf tlX \_. V Jf VJJ » ij i V \^ III 

l Juneauetal. 1 |2005t | Le Floc'h et alj |2005D and various signs of observation time if spectroscopic redshifts are used to derive the 

this downsizmg are now seen. But precise measurements of the luminosities of the sources. Photometiic redshifts can also be 

rate of stellar formation occurring at high redshift are still needed ^^^^ complement spectroscopic redshifts, but their uncertain- 

to efficiently challenge the latest models of galaxy formation. In ^-^^ well-quantified only at low redshift, where specti-oscopic 

other words, additional constraints on the modehng of the evolu- j-edshifts are available to calibrate them. The second reason is 

tion of the cosmic SF history should be inspired by observations, ^j^^^ k-corrections of 24 /zm light becomes large and hazardous 

Recently, very deep surveys were designed to probe SF in at z > 2 where the restframe wavelength falls in the PAH fea- 

the distant universe. For instance, mid-infrai-ed hght (at 15 and ^^j^es. The derivations of total infrared luminosities and SFRs 

24fim) collected by the ISO and Spitzer telescopes has been ^re therefore uncertain. In the following, we call this approach 

"direct method" and it consists in deriving luminosity functions 

Send offprint requests to: Damien Le Borgne, e-mail: and SFRD from mid-IR light collected in redshift surveys. 
leborgne@iap.fr 
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However, these studies present some severe limitations. 
First, a full multi-wavelength approach has not yet been used 
to measure these quantities. Indeed, the studies cited above only 
use observations in one mid-IR band to extrapolate to a total IR 
luminosity and derive an SFR from uncertain calibrations (e.g. 
iKennicutl 1 1 9981) . Moreover, various depths and areas must be 
explored simultaneously. On the one hand, only very deep sur- 
veys are able to probe the lowest levels of SF in distant galax- 
ies, which is necessary to account for the total volume-average 
SF activity potentially dominated by numerous galaxies with 
low SFRs. But these very deep surveys necessarily probe only 
a small area in the sky. On the other hand, large (and therefore 
shallow) surveys are also needed to probe the populations of 
sources presenting a low density on the sky. It is the case for 
very low-redshift sources and it might also be the case for the 
distant populations of ULIRGSQ for instance. 

A natural way of exploiting this multi-wavelength and 
multi-scale information all-together is to adopt a global mod- 
eling approach. Ideally, the models (defined as the combina- 
tion of a library of IR SEDs and of an evolving LF) must 
be able to simultaneously account for all the counts observed 
at all IR wavelengths, from faint to bright sources. In addi- 
tion, they must also account for the constraints brought by 
measurements of the CIRB, and from LFs mea s ured w ith the 
"direct method". For instance, ICharv & ElbazI (|2001') (here- 
after CEOl) and [Lagache et al.1 l2003h (hereafter LDP03) or 
iFranceschini et al. (l2001h have found models that are able to re- 
produce most of these constraints. However, some adjustments 
of the LFs or even of the SEDs by these authors were needed 
to reproduce the most up-to-date observations. This modeling 
approach is powerful but it is also subject to caveats. Indeed, 
some important choices must be made for how to parameter- 
ize the shape (e.g. a double power law, a Schechter function, or 
the local 15 yum LF from Xu (2000) converted into Lir) and the 
evolution of the LFs (e.g. with factors (1 + z)"'- and (1 + z)"")- 
These parameterizations rely mainly on physical intuition and 
sometimes require adding more degrees of freedom. This is the 
case for the two populations of sources introduced by LDP03, 
with their LFs evolving separately from each other Moreover, 
these models cannot claim to be the only possible representa- 
tion of the true cosmic LFs or of the true SEDs. They are gener- 
ally good enough to reproduce current observations, but they are 
never provided with a range of uncertainties. 

In this paper, we add new constraints to the decomposi- 
tion of the "Lilly-Madau" SFRD diagram using a powerful non- 
parametric inversion technique that blindy and simultaneously 
exploits the information from the published multi-wavelength 
IR galaxy counts in deep and small, as well as large and shal- 
low, surveys. From these counts alone, and without any input 
information on redshifts, we derive the range of all possible evo- 
lutions and shapes of the IR LF. 

In Sect. 121 we present the inversion method in the general 
case (see also Appendix|A]i. We present in Section|3]the data and 
our choice for an SED library used in this study. Section |4] con- 
tains our results: the counts inversion and the corresponding IR 
LFs, together with the inferred cosmic SF history. In Section |5j 
we validate our inversion by comparing our empirical modeling 



' We adopt the following notation in this paper: 
"Normal" galaxies have Lir < 10" Lq 
"LIRGS": 10"Lo < Ljr < IO'^Lq 



of counts, LFs, SFRD, and stellar-mass density evolution to bib- 
liographic data (see also Appendix |B] for the robustness of the 
inversion). Finally, we discuss the results and give our conclu- 
sions in Section |6] Some predictions for forthcoming Herschel 
observations are also given in AppendixICl We use a cosmology 
defined by Ho= 70 km s" ' Mpc ', A= 0.7, ^^=0.3. The IMF 
is assumed to be lSalpete^ (fT955l) unless otherwise stated (i.e. in 
Sect.El. 



2. Method: non-parametric inversion of deep multi-/! 
IR galaxy counts 

In this section, we present a new approach to infer the evolution 
of IR LFs of galaxies from multi-wavelength and multi-scale ob- 
servations of galaxy counts. It consists in a phenomenological 
modeling approach, similar to the works of CEOl, LDP03, and 
otherfl However, the modeling is made here in a global and 
flexible way: assuming that a given library of SEDs is able to ac- 
count for the spectra of galaxies at any redshift (an assumption 
tested from the limits of the method's success), we search blindly 
for all shapes and evolutions of the total IR LF that are able to 
reproduce the multi- wavelength IR counts and the CIRB. This 
method is non-parametric; i.e., it does not depend on a param- 
eterization of the LF (in shape or evolution). It must be noted, 
however, that the underlying model of SEDs (on which the in- 
version depends) can, itself, involve one or several parameters 
that must be fixed for the inversion. 

Our method exploits data from infrared surveys designed to 
probe high-redshift populations by using their observed galaxy 
number counts. However, the redshifts of the sources are not re- 
quked, which makes this method quite versatile0. 

An important advantage of this technique, one that makes it 
different from all previous models for IR galaxy number counts, 
is that it provides an automatic, hence objective, way of sam- 
pling the range of possible histories of the IR luminosity func- 
tion. While previous studies have always presented their favorite 
model for fitting the IR galaxy counts, the present work spans 
the range of all possible evolutions that are consistent with the 
observations. This results in two major improvements. First, the 
modeled cosmic SFR history or luminosity functions per redshift 
bin are presented with their error bars. Second, it will allow us to 
discuss the limitation of local IR SEDs at reproducing the prop- 
erties of distant galaxies. Indeed, if after spanning all possibili- 
ties, we still find that the fit is not complete, this will demonstrate 
that the IR SEDs must be revised, either because they provide an 
incomplete description of local galaxies or because they evolve 
with redshift. 

In this section, as well as in appendix A, we present the tech- 
nical and mathematical aspects of our counts modeling. 

2.1. Linear matrix modeling of fhe counts 

The deep galaxy counts can be seen as the projection on a 
flux scale of the SEDs of galaxies of various luminosities. 



"ULIRGS' 
"HLIRGS' 



L,R > lO'^Lo 



- This approach is at variance with the classical "direct" methods 
used for this purpose, in which luminosities of individual sources are 
derived from their measured redshifts and fluxes in a single mid-IR 
band. 

^ It is worth noting that we can take (optionally and when available) 
the redshift information into account by using as priors the luminosity 
functions measured at low redshifts from direct methods. We will show 
in Sect.|4]that this knowledge of the luminosity-redshift distribution at 
low redshift (z < 2) brings actually little new to our results. 
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masses, types, etc., distributed in redshift. Therefore, to repro- 
duce counts, one needs at least a description of the number of 
galaxies at various redshifts per unit volume, the SEDs of these 
galaxiefl and a cosmology. 

In the following, to describe the numbers of galaxies in space 
and time, we use total infrared luminosity functions (LFir). 
Doing so, we simultaneously assume that the SED of a galaxy 
seen in the IR can be efficiently described by the sheer knowl- 
edge of its redshift and Lir. Although such a description can 
be regarded as simple-minded, the current knowledge of galaxy 
SEDs in the IR is not much better than an empirical description 
parameterized only by Lir. Therefore, the model description that 
we have chosen is suited to the current limitations of our under- 
standing. 

For a galaxy (which IR SED is known) lying at redshift zo 
with a luminosity LjJ^, we can easily compute a k-correction and 
a distance modulus to obtain the flux density 5o that one would 
measure in a given IR filter centered at the wavelength A^. This 
conversion only depends on the cosmology and on the SED of 
the galaxy. 

A luminosity function being only the description of the num- 
ber of sources per comoving Mpc^ as a function of redshift and 
infrared luminosity Lir, we can then compute a matrix that con- 
verts this evolving LFir into numbers of galaxies seen at various 
fluxes and through different IR filters. This matrix is simply the 
response function of the conversion from (z, Lir) to {S,A) pre- 
sented above. Again, this matrix depends only on the cosmology 
and on the library of SEDs. Since it includes k-correction, dis- 
tance eff'ects (dimming as a function of the square of the lumi- 
nosity distance), and redshift effects (flux density stretching and 
dimming), we call it the "k-i-d" matrix in the following 

After discretization, and using a matrix notation, the inverse 
problem can be formalized as 

Y(^, S ) = M(^, 5 , z, Lir) ■ X(z, Lir) , ( 1 ) 

where Y is the matrix containing the number counts at fluxes 
S in bands A, M is the above-mentioned "k+d" transformation 
matrix, and X is the LFir that is a function of only z and Lir (see 
Appendix A for details). 

Therefore, our problem involves inverting this linear equa- 
tion to find the evolving LFs (i.e. redshift dependent number 
counts per unit volume), X, from the known values of Y (the 
wavelength dependent observed number counts per unit area). 
Because the matrix, M, is not square, and because the num- 
ber counts are noisy and must be positive, the solution is not 
quite as simple as using the pseudo inverse: X = M*"'' ■ Y, and 
requires computing a regularized solution as discussed bellow 
and explained in Ap pe ndix A. The reader ma y also refer to, e.g. 
iPichon et all (120021) or lOcvirk et"aLl (l2006 b'). The uncertainties 
on the observed counts Y are taken into account through an ad- 
ditional error matrix W that makes it possible to compute a 
between the model X and the data Y and to derive uncertainties 
on the recovered LF. 

In this formalism, we choose to describe the luminosity func- 
tion in a logarithmic scale, taking logjoLiR instead of Lir every- 
where. This is justified by the better conditioning of the inver- 
sion in this case because the LFir generally spans several orders 
of magnitudes in luminosity. For the same reason, the counts 
are treated numerically through their Euclidian differential form 

In this subsection presenting the formaUsm in the general case, an 
evolution of the SEDs is considered as possible. 

^ If the SEDs vary not only with Lir but also with z, we could call it 
the "e-l-k-Fd" matrix. 



dN/dS X S^'^ (units of mJy'^ deg"^) which varies slowly with 
flux. The evolution is measured as a function of logio(l + z). 
Indeed, in the following, we impose a smooth evolution in red- 
shift of the luminosity function and we need to define the tempo- 
ral parameter on which this smoothing applies. We find that the 
simplest redshift description that corresponds roughly to a regu- 
lar time sampling is actually logj,)(l -i-z). A parameterization with 
z or with log 10 z would leave too much room for strong variations 
in LFir at early and late times, respectively. More quantitatively, 
we discretize the problem in bins regular in (51ogiQLiR=0.1 and 
in 51ogio(l +z) - 0.015, which corresponds to the sampling that 
is good enough to produce counts with regular flux sampling of 
51ogioS-0.1. 

2.2. Regularizing the inverse probiem 

In practice, provided the number of bins in LFir is large enough, 
this problem is "ill-posed": there is possibly a large number of 
LFs that are able to satisfy Eq. ([T]i perfectly, hence overfilling the 
noisy counts. However, many of these solutions are unlikely and 
not physical so that we need to "regularize" the problem to ob- 
tain valid solutions. The first natural constraint is the positivity of 
the LFir (there are no such things as negative numbers of galax- 
ies). This imposes an iterative, CPU-costly approach to the prob- 
lem, and a choice for an initial guess (see below). Moreover, we 
want to avoid solutions that are not meaningful given our noisy 
finite set of counts, such as LFs that are chaotically varying as a 
function of Lir or redshift. Therefore, we penalize the inversion 
with an extra term, added to the formal x^, which enforces the 
smoothness of the L Fir both in Lir and in z (see Appendix A and 
lOcvirk et al.ll2006b, for the details of the formalism). With these 
constraints, the solutions, X (which depend on the choice of the 
initial guess required for the nonlinear optimization), are reason- 
able and can serve as a solid basis for future works. We can also 
optionally impose external constraints as supplementary priors, 
namely the low-redshift luminosity functions obtained from di- 
rect methods. 

To obtain the range of all realistic solutions for the nonlin- 
ear optimization problem, we explore a wide range of random 
initial guesses, in a Monte-Carlo approach with at least 100 real- 
izations. The range of LFs spanned by the initial guesses is much 
wider than the range of the final LFir (X) obtained after conver- 
gence, which lends credibility to our study's completeness 

Finally, some of the solutions do not match the constraints 
brought by CIRB measurements. We filter out these invalid so- 
lutions, a posteriori, leaving only the most realistic luminosity 
functions. 

3. Data and SEDs used in this work 

3.1. Input multi-wavelength counts 

The multi-wavelength counts that we invert in this work are pre- 
sented in Fig.[T] They summarize bibliographic observed galaxy 
counts at 15, 24, 70, 160, and 850//m. This compilation is not 

* We could also have used a more classical approach to character- 
ize the uncertainties in the recovered LFs by computing the posterior 
variance covariance of the parameters away from the itterative nonlin- 
ear solution, and compared it to the corresponding initial prior (hence 
constructing the so called information matrix). These uncertainties nat- 
urally arise from the uncertainties on the observed counts through the 
error matrix W. Instead, we favor here Monte-Carlo simulations be- 
cause they yield similar - if not more robust (since we span ranges of 
possible nonlinear solutions) - constraints on the recovered LF. 



4 



Le Borgne et al.: Cosmic star-formation history from galaxy counts 



LW3/IRS 15/16/i,m 



o 1000: 



c/i 
* 



-3' 



Elboz et ol. (1999) 
I Gruppionl et ol. (2002) 
Q from Teplitz et ol. (2006) dotal 



100 : 



1.0 

Flux (mJy) 



10.0 



100.0 



1000 



100 



MIPS 2Vm 



Chory et ol. (2004) 
I Popovicti et ol. (2004) 
Q Le Floch et ol. (2005) 



0.1 



1.0 
Flux (mJy) 



10.0 



100.0 



MIPS 70yLim 



10° 



"tn 10': 



10° 



A Froyer et ol. (2005o) 

I Prayer el ol. (2006b) 

Q Dole et ol. (2004) 

m Dole stocking (priv. comm.) 



-r ? 5 %$T 



^1 ^ 



10 

Flux (mJy) 



100 



^ 10000: 

J, 

""^ lOOOt 
(/) 

100 . 



SCUBA 850yLtm A Coppin et ol. (2006) 



I 



Flux (mJy) 



10° - 



10- 



10" 



ISO/MIPS 170/160;tim A D°ie et oi. (2001) 

I Dole et ol. (2004): Morano 

Q Dole et ol, (2004): CDFS 

^ Dole stocking (priv. comm.) 

T 5 



10 100 
Flux (mJy) 



1000 



Fig. 1. Observed 15/16jum, 24/vm, 70yum, 160 /im, and 850//m counts of galaxies from various surveys. The shaded area rep- 
resenting 1 cr uncertainty area are determined semi-automatically from the data points (smoothed envelopes of the data points 
after excluding the str ongest out liers, particularly at 15 /im), and are used as an input for the inversions. Instr uments/fields: ISO 
LW3/various fields (Elbazetal. 1999), ISO LW3/ELA1S-S (Gruppioni etaTI 120021). Spitzer IRS/GOODS ( ITeplitz et al.ll2006l) . 
Spitzer MlPS24/GOq DS tCharv et al., .2004.) , S pitzer MlPS24/vario us fields dPapovich et al., .2004) , Spitzer MIPS24/G OODS 
(iLe Floc'h et al.l l2005h . S pitzer MIPS70/F LS (iFraver et al.l l2006ah . Spitzer MIPS70/GOODS-N (iFraver et al.1 l2006bh . ISO 
ISOPHOT170 /FIRBACK jPole et al.1 l200lb , Spitzer MIPS160/Marano-i-CDFS jDole et alj l2004l) , JCMT SCUBA/SXDF-hLH 
dCoppin et al.ii2006.) . 



exhaustive but we retained significant surveys for each IR band, 
mixing wide and shallow ones to avoid cosmic variance effects 
as much as possible, and deep ones for faint source detections. 
The data shown in Fig. [T] are actual observations, except for the 
70 and 1 60 Aim points fro m H. Dole (private communication, 
and also iDole etalJl2006l) . which are obtained from stacking 
24//m sources in several flux bins and using an MIR-FIR ob- 
served correlation. Although these points are not strict measure- 
ments, the MIR-FIR relationship is tight enough (down to less 
than 0.1 mJy) to be confident in the stacking results. We group 
the 15 and 16/im counts, as well as the 160/1 70//m counts, thus 
neglecting the small differences to the counts (a few percent) that 
are caused by the slightly different k-corrections. Our inversion 
is applied to this compilation. 



3.2. Library of SEDs 

The IR SEDs that we use in this study are taken from the empir- 
ical library of CEOl, which defines a bijection between Ljr and 
the IR SED from 3 to 1000 yum. Although this library is based 
on correlation observed in local galaxies, we suppose in the fol- 
lowing that the SEDs of galaxies with a given Lir do not change 
with redshift. It does not necessarily mean that SEDs of indi- 
vidual galaxies do not evolve, but that they must evolve along 
the local SED-Lir correlation. It is worth noting that when a de- 
scription of the evolving IR SEDs of galaxies becomes available 
(e.g. with Herschel), the technique that we use can be applied to 
these evolving SED libraries. 

Moreover, we 'clip' this library for the highest IR luminosi- 
ties, using the Lir=10'^^ L© SED shape for galaxies more lu- 
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4.1. Counts inversion 
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Fig. 2. 24 yum vs total IR luminosities (or the SFR equivalent of 
these luminosities assuming a Salpet er IMF). The shaded gray 
area represents the observations of Papo vich et alj (2X)01) Fig. 
11, considering that most of this intense SF is probed by the 
IR light, neglecting the UV contribution. The solid lines corre- 
spond to the "clipped" CEOl library for a set of virtual galaxies 
at redshifts z=1.5, 1.8, 2.2 and 2.5 (from left to r ight) with fluxes 
fainter than 1 mJy, for which the Pap ovich et al.l ( 1200 7) observa- 
tions are valid. A pure CEOl library would have followed the 
one-to-one slope. 



minous than this value. This makes the SEDs colder than in 
CEOl for ULIRGs and makes it compatible with most of the 
known data, from low-redshift, low-luminosity to hi gh-redshift 
luminous galaxies (see CEOl. IPapovich et al.l l2007l) . Figure |2] 
presents the relation between 24fim and total IR luminosities 
that is predicted b y the "cli pped" CEOl library, together with ob- 
servations from dPapovich et alj|2007l) in the range z = 1 .5 - 2.5. 
This clipping can be justified by the SEDs of high luminosity 
galaxies (ULIRGs) like those seen at high redshift being poorly 
known in the local universe: their dust temperature can only be 
measured by Herschel. Therefore, the CEOl SEDs are extrapo- 
lated within this luminosity range. We chose to clip the IR dust 
temperatures of the CEOl library to those of the luminosity range 
really observed in the local universe, rather than extrapolating 
them. 



4. Results 

In this section, we present the counts, luminosity functions, and 
SF history that are modeled from the counts inversion, as well 
as the effect of using or not using low-z luminosity functions as 
priors. We compare these results and their redshift decomposi- 
tions to measurements obtained from direct methods and biblio- 
graphic data only in Sect. |5] 



This inversion model has been designed to reproduce infrared 
galaxy counts and as a result to derive a range of possible to- 
tal IR LF as a function of redshift, that can be converted after- 
wards into a range of cosmic SF histories. The success of the 
model can be visually tested by comparing the range of predicted 
galaxy counts with the observed number counts and their disper- 
sion (see Fig.|3]l. 

At first glance, one can see that the observed counts are well- 
fitted over the whole IR range, from 15 to 850 jum. In particu- 
lar, the bumps at 15 and 24yum are reproduced simultaneously. 
The 850yum differential counts are well-fitted too: negative k- 
corrections make it possible to see a high-redshift population of 
galaxies, namely ULIRGS and HLIRGS (Lir> IO'^Lq) at z > 2, 
which are hardly seen at other wavelengths except in faint 24 
and 70 /im counts. 

Interestingly, one can see from Fig. [3] that, while the only 
strong constraint that is imposed on the model is to keep a 
smooth dependence of the LF with redshift and luminosity, the 
model is unable to perfectly fit the observed number counts and 
their dispersion at all flux densities and wavelengths even though 
the whole range of possible LF and associated redshift evolu- 
tion has been spanned blindly. Some solutions tend not to fit 
the 15fj.m counts perfectly, the 70jum counts are slightly over- 
produced, and more strikingly, the 160/im counts are underpro- 
duced around 20 mJy. 

A major strength of this model is to provide an objective and 
statistically significant way to test a given library of template 
SED. Indeed, the discrepancy between the model and observed 
counts cannot arise from the LF itself since it was allowed to 
vary both wit h luminosity and reds hift with a high degree of free- 
dom (see also lOcvirk et al.ll2006a . for a discussion on the corre- 
sponding biases in a slightly different context). It must therefore 
arise from the library of template SEDs that is used as an in- 
put for fitting the number counts (through the "k-i-d" matrix M). 
The CEOl Ubrary of template SEDs that is used here reflects 
the median trend of local galaxies and was found to be statisti- 
cally consistent w i th the radio-infrared correlation up to z ~1.3 
(lElbaz et al.lll"999t lAppleton et al]l200 4). with the mid-infrared 
observations of galax ies up to z ~1 with moderate variations 
(iMarciUac et al.ll2006|l. and with massive galaxies selected with 
the BzK technique (iDaddi et al.l l2007b). The origin of the dis- 
crepancy can therefore come from three possible origins. The 
first two possibilities compatible with no evolution of the in- 
frared SED of galaxies are (i) a bias towards cold galaxies due to 
the shallow depth at 160yL(m within the dispersion already exist- 
ing at z ~0; (ii) the existence of a subpopulation of cold galaxies, 
akeady present locally but not yet identified due to limited con- 
straints on both sides of the peak emission in the far infrared. 
A third possibility would be that the i nfrared SED of galaxi es 
evolve as a function of redshift (see e.g. lChapman et all 20021) . 

It is not possible to disentangle between the three possibili- 
ties based on the present dataset. However, a forthcoming paper 
will study this issue in detail using a stacking analysis at 160 //m 
(Magnelli et al., in preparation). It must be noted though that 
the method described here is versatile enough to allow its user 
to test any library of template SED against existing constraints 
from galaxy counts and the infrared background. It will there- 
fore be a straightforward matter to check whether any change or 
evolution in the SEDs can reproduce the number counts at all 
wavelengths and flux densities. We refrained from making these 
adjustments to the counts at this stage since any of the previously 
mentioned alternatives is equally possible. 
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Fig. 3. Mid- and Far-IR counts obtained from the non-parametric inversion of the observed counts. The data counts that are the basis 
of the inversion are represented by the shaded blue zones (fitted values ± 1 cr, corresponding to the zones of Fig.[T]i. The best fit to 
the counts is shown by the solid red lines, together with the range of allowed solutions (red dashed areas). The vertical gray zones 
delimit the flux ranges where some counts are observed. All the counts modeled here satisfy the CIRB constraints. The oblique 
dashed line corresponds to a statistical hmit of a one square degree survey (less than 2 sources per decade in flux). 



4.2. Solutions: Range of evolving luminosity functions 

Looking now at the main output of the model, i.e. the evolution 
of the total IR luminosity function with redshift, we note that 
large parts of the LF are constrained very well by the inversion 
(see Fig. |4]i. In particular, the number density of galaxies with 
10.5 < logio^iR < 11.5 at z < 0.5 is tightly constrained, as 
are the numbers of galaxies with higher luminosities at higher 
redshift. 

We will present in Sect. 15. 21 a comparison of these solutions 
to the LFirS obtained from direct measurements at low redshift. 

But before doing so, we note that a fraction of the solutions 
present a knee in the LF, particularly at 0.3 < z < 1 around 
10'" - 1O"L0, which is partly responsible for the bump seen 
in the number counts at 24 and 15 yum at the corresponding flux 
densities. However, not all solutions of the inversion technique 
present such a strong knee, which might be seen as an artifact. 

4.3. Evolution of the star-formation activity 

The SFRD can now be estimated from the LFs that were ob- 
tained from the counts inversion. The total infrared LFs are inte- 
grated over the whole range of luminosities down to 10^ L©. The 



resulting total infrared lumi nosity density i s then converted into 
a SFRD using formula ^ (IKennicutllll998h : 

SFR [Mo yr-'] = 1.72 x lO"'" Lir [Lo] . (2) 

Therefore, the regions shown in Fig. Ineffectively represent the 
range of all possible SF histories that are compatible with the 
multi-zl counts and the CIRB. The SFRD history obtained from 
the LFiR that produces the best fit to the counts is also shown for 
each luminosity class of galaxies. 

This inversion shows that, indeed, an IR downsizing is at 
work: "normal" galaxies dominate the SFRD at low redshift (al- 
though the contribution of ULIRGS and HLIRGS are poorly 
constrained in the low-z range because they would correspond 
to bright and very rare sources, not easily seen in the counts of 
current deep surveys). At z > 0.8, LIRGS dominate the SFRD, 
whereas the contribution of ULIRGS peaks at z > 2. These re- 
sults will be compared in detail to measurements from direct 
methods and bibliographic data in Sect. 15. 31 

As mentioned before, it is also possible to use an additional 
prior for the inversion: direct measurements of the 8 or 15 jum LF 
below z = 2. Although subject to many caveats (e.g. the strong 
dependency on the PAH modeling of the SEDs), they can be 
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Fig. 4. All possible solutions for the evolving LFs that best reproduce all the IR counts of Fig. [3] and the CIRB constraints. As in 
Fig. [3] the thick red line and the dashed red areas represent the best-fitting solution and the range of allowed solutions, respectively. 
By construction, no conclusion on the LF can be made in the vertical gray-shaded areas where a divergence is expected because 
these objects are not seen in the counts because of the flux limits of current IR surveys. At the top of each panel, the ranges of 
LiR probed by IR counts are shown as horizontal lines, from 15 to SSOyKm (bottom to top). Diamonds, triangles, and filled circles 
correspond to 0.1, 0.3, and 3 mJy fluxes, respectively. Bottom right panel: range of CIRB corresponding to the inverted LFs. Data 
points (in gray) are taken from the compilation bv .Dole et al.iC2006.) . 



used as a prior to guide the inversion, and at least constrain the 
solutions at low redshiffl 



^ Technically, this supplementary prior involves again penalizing the 
formal by adding an extra term that measures the distance be- 
tween the observed LF and the solution (this corresponds to /j 9^ 
in Appendix [A]). The LFs that are too different from the known low- 
redshift LF are therefore strongly penalized, so excluded de facto. 



We checked the eff^ect of using this prior from direct mea- 
surements at low redshift. We observe that some uncertainties in 
the LIRGS and ULIRGS contributions are slightly tightened, but 
we also note that the trends are basically unchanged. We inter- 
pret this surprising result as follows: the leverage that we have 
access to by inverting galaxy counts on a very wide wavelength- 
basis (from 15 to 850yum) is large enough to provide a realis- 
tic description of the redshift distribution of the sources on a 
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Fig. 5. Star-formation rate density since z = 5 derived from the counts inversions. The SFRD is obtained from the range of all 
possible luminosity functions derived from the counts and respecting the CIRB constraints. The solid lines correspond to the best 
fit to the counts and the transparent shaded areas show the range of uncertainties. From top to bottom at z=0.8: blackball galaxies, 
green=LIRGS, blue=normal star-forming galaxies, orange=ULIRGS, red=HLIRGS. 



statistical basis. This is likely to only be possible because the 
library of SEDs that we use seems close to the real SEDs of 
galaxies, on average (again, in the statistical sense) at any red- 
shift lower than z - 3, thus avoiding a complete blurring of the 
de-projection of multi-wavelength galaxy counts onto the lumi- 
nosity function space. Therefore, various populations of galax- 
ies at diff'erent redshifts are seen at various wavelengths, which 
considerably reduces degeneracies and enables us to recover the 
history of IR galaxies as a whole. 

Since the philosophy of this paper is to remain as conserva- 
tive as possible, we choose not to use the low-redshift LF mea- 
sured from direct methods as a prior in the following. Indeed, 
doing so would only slightly change our results, and it would in- 
troduce a source of potential additional errors propagating from 
the errors intrinsic to direct methods (k-corrections in the mid- 
infrared or redshifts). 



5. Validation: comparison with direct 
measurements at low redshift 

After presenting the global outputs of the model in the previ- 
ous section, we detail here the redshift decomposition of the in- 
version results and we compare them to measurements obtained 
from direct (redshift-based) methods. This comparison is partic- 



ularly challenging since no redshift information was used as an 
input in the inversion. 



5.1. Redshift decomposition of tine mid-IR counts 

To validate our inversion results, we need to compare the red- 
shift decomposition of the IR counts to observations from a 
direct method. To do so, we compared our results to data 
from the GOODS survey (RI. Dickinson for GOODS-Spitzer, 
RI. M.Giavalisco for GOODS-HST) originally presented in 
iGiavalisco et al.l (l2004t) . This survey consists of two fields which 
have been subject to several studies at various wavelengths in 
the past few years. We investigated the redshift decomposition 
of the counts at 16, 24, and 70yum using the optical counter- 
parts of these sources in GOODS and making use of the spectro- 
scopic and photometric redshifts. The sample that we used cov- 
ers a total area of 0.07 square degree on the sky. Although this 
area is quite small and cosmic variance might affect our study, 
we found that the luminosity functions measured from the di- 
rect method are similar in both GOODS fields, making them 
compatible within 2cr. The spectroscopic completeness is high 
(60% at S 24^/,, - 30 yuJy for z < 1 .5 sources), and we comple- 
mented tJiem_wA_£hoto with the code 
Z-Peg (iLe Borgne & Rocca-Volmerangell2002h with a precision 
Az/(1 +z) - 0.04 to 0. 1 depending on the redshift of the sources. 
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Fig. 6. Top and bottom left: Counts at 16, 24, and 70jum modeled (solid lines) from the inversion, on the one hand, and observed 
(dotted lines) in GOODS fields, on the other. The counts are decomposed in redshift bins (black=all redshifts, blue=z < 0.5, 
green=0.5 < z < 1.5, orange=1.5 < z < 2.5, red=z > 2.5). Bottom right: CIRB decomposed in redshift (same color-coding). The 
oblique dashed line corresponds to the limit in statistics due to the smallness of a field like GOODS North-t-South or 0.07 square 
degrees: less than 2 galaxies per flux bin of width 51ogF=0.1 dex are expected below this limit. 



We then used the Vmax formalism to correct from incomplete- 
ness at the lower flux limit. The galaxies showing signs of AGNs 
(identified from X-rays or optical emission lines), were excluded 
from the sample of 24 /im sources. Doing so enabled us to use 
SEDs of galaxies to compute k-corrections and only slightly af- 
fects our results, mainly at the very high-luminosity end of the 
luminosity function at moderate (z = 1) or high (z > 2) redshifts. 

Figure |6] presents the resulting comparison of the redshift 
decomposition of the counts obtained from direct and inverse 
methods. At 24 /im, our best solution for the recovered LFir 
indeed produces a redshift decomposition of the counts that is 
compatible with the observed ones. The match of 15 fim counts 
as a function of redshift is poorer because CEOl templates rep- 
resent the fl uxes at this wavelengt h less well for 0.5 < z < 1 
galaxies (e.g. lMarcillac et al.l2006h . As for the 70 yum counts, the 
observed decomposition is not complete at faint fluxes, making 
the comparison hazardous. 

5.2. Comparison to direct measurements of tiie infrared 
iuminosity functions 

We now compare the range of luminosity functions LFir ob- 
tained from our inversion to some measurements of the LFir ob- 
tained from the direct method. As noted before, several studies 
(iLe Floc'h et al.ll200l iBabbedge et aTll2006l: ICaputi et al.ll2007l) 



have measured the 15 or 8//m luminosity functions, which can 
be converted to LFir if a library of SEDs is ass umed. We show in 
Fig.|2]the results from jLe Floc'h et alj|2b0 5') for reference. We 
also provide on the same figure our own direct measurements 
of the evolving LFir that we derived from 24 /im galaxies seen 
in the GOODS fields (North-i-South). We chose to do these di- 
rect measurements again for two reasons. First, for consistency 
with our inversion, we used the same library of templates for the 
k-corrections (used for MIR to total LIR conversion^ Second, 
our data at 24 fj.m reaches a depth of 24 /iJy instead of 80 /iJy in 
previous works, extending the LF to the faint end. The complete- 
ness at this depth is 85% (Charv 2007). Not surprisingly, our di- 
rect measurements of the LFir from 24 jum fluxes are compatible 
with previous studies that used the same niethod . In particular, 
we find the same values as lLe Floc'h et al.l (l2005l) who explored 
the z < 1 domain. However, we reach higher redshift galaxies 
thanks to the depth of our sample. However, k-corrections are 
very uncertain for these high-z sources (24 /im corresponds to 
less than 8 fim at z > 2, a range where the library of SEDs is not 



^ Note that lLe Floc'h et il\ d2005l) errors bars are rather large because 
they partly reflect the uncertainties in the SEDs of galaxies. We checked 
that using another hbrary of SEDs (e.g. the LDP03 h brary) produces 
differences in our LFs comparable to the error bars of iLe Floc'h et al.l 
1 12005) . 
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Fig. 7. Total IR luminosity functions derived independently from the counts inversion and from direct measure ments. "Clipped" 
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validated). Therefore, we use these high-redshift measurements 
with caution, as mentioned earlier 

It is striking in Fig. |2]that the LF obtained from the direct 
method is consistent with the best-fitting LF derived from the 
counts inversion in the common range that they probe. This is re- 
markable because no information on the redshift of the sources 
was used in the inversion. This means that all the constraints 
that one can get from these direct measurements are not really 
needed for the inversion: although we do not use the redshift of 
the sources, we recover the observed redshift distributions that 
are here expressed equivalently in terms of evolving luminosity 



functions. The interpretation is the same as for the redshift de- 
composition of the counts, which are different views of the same 
phenomenon. 

The meaningfulness of the agreement between the observed 
and the recovered LFs is strengthened by the fact that that both 
Lir luminosity functions are measured or estimated using the 
same library of SEDs, thus using the same k-correction in a con- 
sistent way. Using another library of SEDs for both methods 
would produce the same kind of agreement, although the pre- 
cise shape of the LFs would be slightly different from what we 
obtain here with the "clipped" CEOl library. 
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5.3. Cosmic star-formation history vs literature 

In this section, we compare the SFRD history that we derived 
from the multi-/l counts inversion to "direct" measurements at 
low redshift. 



5.3.1. Total SFRD 

The comparison of the total SFRD derived with this inversion 
tech nique ma tches the corn pilation of direct measurements from 
iHop kins & B eacomI (|2006|) (Fig.[8]l. Moreover, we also show in 
Fig.[8]the SFRD that we derived from 24 ^/m data in the GOODS 
fields. We must stress here that these measurements, obtained 
with the "direct" method, are independent of the SFRD inferred 
by the counts inversiorQ. We only present them here for com- 
parison to the inversion results in this subsection. The total in- 
frared luminosities measured for individual 24 fj.m sources in the 
GOODS fields were summed up in redshift bins and converted 
to SFRDs. Like severa l authors before us (e.g.lFlores et al.ll999t 
iLe Floc'h et"an 120051: iPerez-Gonzalez et all l2005l) . we could 
have extrapolated the measured luminosity functions in the faint 
end to obtain the total SFRD. But this is somewhat dangerous 
and depends strongly on the parameterization of the LF fits, 
so we chose not to extrapolate the luminosity functions at faint 
fluxes to estimate the SFRD. Instead, we used only the sources 
brighter than our flux limit at 24 fj.m (24 yuJy) to estimate lower- 
limits for the SFRD at every redshift. Therefore, all the points in 
this figure should be considered as lower limits. 

Our inversion technique, on the other hand, makes it possi- 
ble to partially avoid such caveats. First, unlike direct methods, 
more than one band is used to estimate the total IR LF. Second, 
the extrapolation of the LF in the faint-end is achieved automat- 
ically via the only constraint of having a smooth variation in Lir 
and z. Therefore, the area of uncertainties we propose here are 
very likely more robust than previous estimates because they use 
more data and are less dependent on parameterization for both 
the shape of the LF and its evolution. 

The SFRD obtained from the counts inversion and our mea- 
surements with the "direct" method in GOODS are in rather 
good agreement, which tends to give credit to the inversion. 
One might note, however, that the data points are systematically 
lower than the inversion results, which illustrates the choice of 
not extrapolating the measured LFs at faint luminosities to esti- 
mate the SFRD. 



5.3.2. Luminosity decomposition of tine SFRD 

Figure |9] presents a detailed comparison of the history of the 
SFRD inferred from our inversion and decomposed in luminos- 
ity classes to what can be independently obtained from the "di- 
rect" method with 24yum sources in GOODS. 

Let us first consider the reliability of the direct measurements 
used for comparison. Because we chose not to extrapolate the 
LFs at low luminosity, several points in this figure must be con- 
sidered as lower limits. It is not the case, however, for the z < 1.2 
points for LIRGS and the 1 < z < 2.5 points for ULIRGS. 
Overall, we are limited at very low redshift by small statistics 
and at high redshift by flux limits. Moreover, the contribution 
of HLIRGS should be taken with caution. Indeed, by inspecting 
these luminous high-redshift sources, we came to the conclusion 



' The datasets are independent, indeed. But strictly speaking, it is 
actually not the case for the SFRD measurements because the libraries 
of SEDs used for k-corrections are the same in both methods. 



that a large fraction of the photometric redshifts computed for 
these IR-bright sources have poor precision, leading to catas- 
trophic failures for almost half of the HLIRG sources at z > 1.5. 
This poor performance of the photometric redshifts is not too 
surprising for this class of galaxies at such high redshifts since 
the templates used in the fitting procedure have a relatively low 
level of dust, compatible with most galaxies seen in the current 
optical and NIR surveys. These direct measurements present a 
nice picture of the "IR downsizing", where the cosmic SFRD 
was dominated by brighter and brighter galaxies in the IR as we 
go back in time. Our results confirms the IR view of the cos- 
mic SF histo ry that was explored in previous works up to z = 2 
ICaputi et al l (e.g. '2007). At z = 2, we confirm that ULIRGS 
seem to dominate the budget of the SF activity. 

A number of interesting remarks arise from the comparison 
of the "inverted" SFRD to the "direct" measurements. First, it 
is comforting to see that both methods give consistent views of 
the IR downsizing. In both cases, low luminosity galaxies dom- 
inate the SFRD at z < 0.5, and ULIRGS are dominant at z > 2. 
However, a more detailed comparison of both results provides 
interesting clues to what is really seen in current deep surveys. 
One can easily notice that most data points ("direct" measure- 
ments) are at the lower limit of the area allowed by the inver- 
sion. This means that the 24yuJy limited sample of 24 /zm galax- 
ies probe almost all the cosmic SFRD. However, the same re- 
mark also opens up the possibility that up to 50% of the SFRD 
is not yet resolved in sources down to our flux limit, especially 
for low-luminosity galaxies at z > 0.5. 

Another interesting point is that the inversion does not al- 
low the population of HLIRGS to contribute much to the CIRB 
at any redshift. This somewhat contradicts the direct measure- 
ments for the same objects, which tend to indicate an increas- 
ing contribution of these extreme sources at z > 2. But again, 
we must stress that many uncertainties lie in the observations 
of these distant sources (photometric redshifts, validity of the 
SEDs, contribution of AGNs, etc.). Therefore, we must con- 
clude that this population cannot be too numerous to reproduce 
the deep IR counts, including the 850 yum ones, if the "clipped" 
CEOl SEDs are used at any redshift. We checked that allowing 
more HLIRGS at these redshifts, at a level comparable to the 
observed number, overproduces the 850;um counts. This means 
that either these objects do not exist (and indeed, as we men- 
tioned before, about half of these HLIRGS have a wrong photo- 
metric redshift, hence a wrong luminosity), or the SEDs for these 
objects are very different from the templates in the SED library 
that we use. Of course, both reasons may be at work simulta- 
neously. Interestingly, if the original CEOl library is used (both 
for the inversion and for direct measurements), the situation is 
similar: the SFRD of HLIRGS inferred by the inversion of the 
multi-wavelength counts is still smaller (by a factor of 4) than 
the value inferred from direct measurements of 24 //m sources. 
The main difference, in this case, is that the SFRD of HLIRGS 
are roughly a factor of 5 to 10 larger than what is derived with 
the clipped CEOl library. 

Finally, we notice that the measured SFRD for ULIRGS at 
z < 1.5 (see Fig. |9]l is smaller than the range allowed by the 
inversion. This could be explained by cosmic variance due to the 
small area covered by our GOODS sample. 

5.4. Evolution of the stellar mass density 

In this section, we address the question of the consistency of 
the SFR history that we derive from the inversion model with 
the independent observational constraints existing on its inte- 
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Fig. 8. Total SFRD regions compared to the compilation of direct measurements from iHopkins & BeacomI (|2006|) (empty circles 
and error bars only for z > 1.5 for clarity purposes). The outer light-gray regions corresponds to all possible LFs fitting the counts. 
The inner darker region includes 68% of these models. The dashed line corresponds to the LF producing the best-fit to the multi-/} 
IR counts. The solid triangles with error bars show the integrated SFRD that we measured in the GOODS fields in several redshift 
bins (see also Fig.|9]for a luminosity decomposition of these measurements). 



gral, namely the evolution of the comoving density of stars per 
unit comoving volume. After assuming an initial mass func- 
tion (IMF) and computing the mass of stellar remnants after the 
death of massive stars, it is straightforward to compute the total 
amount of stars that must be locked into galaxies as a function 
of redshift, on the basis of the SFR history. In our computation 
of stellar masses, we account for the recycling of stellar material 
into the ISM, for the mass of stellar remnants (which account for 
about 15% of the total stellar mass at z = 0) and for the evolution 
of the metallicity, using the spectral sy nthesis code PEGASE.2 
jFioc & Rocca-Volmerang3ll997l [19991) . 

This is not the first time that such a computation has been 
performed, but we believe that this is an important test that has 
been the source of discussions in the recent past, in particu- 
lar with the claim that both histories - SF and stellar masses - 
were not consistent with the integral of the SF history producing 
more stars than actually observed at any redshift. Our paper now 
proves this claim to be incorrect. 

Before discussing our own computation, we wish to empha- 
size an important point regarding the effect of the choice of 
a particular IMF in this process. Although various IMFs have 
been proposed in the past, including t op-heavy IMFs for star- 
bursti ng galaxies (see e . g. jE lbaz et all ll995L 119921; iRieke et akl 
119931: iLacev et al.ll2008l:lDavg i2008). no definitive evidence has 
been provided yet for a non-universality of the IMF. The main 
difference that is n ow commonly a ccepted with respect to the 
pioneering work of ISalpeteJ (1 19551) is the finding that the slope 



of the IMF changes around 1 M©, in the direction of having a 
lower contribution of low mass stars to the total mass of stars 
formed, or equivalently a larger contribution of stars more mas- 
sive than 1 Mq (see the review by Cha brier 2003). Nonetheless, 
such a change in the IMF almost equivalently affects both the 
conversion factor used to determine the SFR from LFir and the 
mass-to-light ratio used to derive the stellar mass. In the present 
study, the SFR is derived from the total IR luminosity assuming 
the coefficient computed for a Salpeter IMF by Kennicutt (1998, 
see Eq.|2]l. For a different IMF, such as the Baldry & Glazebrook 
(2003, hereafter BG03) one which shows a flattening below 1 
Mq as discussed above, the SFR would be 0.45 lower and the 
M/L ratio would also be reduced by a similar factor (0.6, com- 
puted using the PEGASE.2 code). 

The evolution of the cosmic stellar mass density with red- 
shift that we computed by integrating the SFR history result- 
ing from the inversion model is found to be in good agreement 
with the latest direct measu rements of galaxy masses of e.g. 
iPerez-Go nzalez et al.l (l2008h . Their published stellar mass den- 
sities were mul t iplied by a factor 0.61 since they were estimated 
with a lSalpeterl (Il955h IMF, which corresponds to the difference 
in mass-to-light ratios in the K band. It is clear that, at all red- 
shifts probed, the range of SF histories that result from the in- 
version technique are consistent with the measured stellar mass 
density. This consistency indirectly reinforces the likelihood that 
the inversion technique spans a reasonable range of possible his- 
tories. 
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Fig. 9. History of the SFRD decomposed in four infrared luminosity classes. The inversion results (full set of models and 68% 
inner region) are shown with gray-shaded areas. The total SFRD history (68% inner region) is also shown for reference in each 
panel as dotted lines. Our direct (and independent) measurements from the GOODS survey are shown as triangles. Empty triangles 
are used for bins affected by completeness. The z = 3 point for HLIRGS is subject to caveats (see text for details) so is probably 
overestimated. 



Finally, we note that after submission of the present paper, an 
erratum was published by Hopkins & Beacom (2006) in which 
they recognize that their computation was erroneous and that, 
contrary to their initial claim, the two histories do not exhibit 
any inconsistency, apart from possibly at the largest redshifts. 
Hence their study is now consistent with our finding, which is 
not surprising since our SFR history globally agrees with their 
compilation (see Fig.[8]l. 



6. Discussion and conclusions 

This paper presents measurements of the evolving infrared lumi- 
nosity function and of the corresponding cosmic SFRD using a 



non-parametric inversion of the galaxy counts in the mid and far 
infrared. 

For the first time, we have exhaustively derived the range 
of possible evolutions for these quantities with a non-parametric 
inversion technique. The input data that were used simultane- 
ously to derive these set of allowed models cover a wide range 
of wavelengths: deep infrared counts observed at various wave- 
lengths (from 15 to 850/im), the cosmic infrared background 
measurements, and optionally the low-redshift measurements 
of the IR luminosity function that are derived from the 24yum 
fluxes. We derived from this inversion the allowed range of SF 
histories, together with the range of stellar mass density evolu- 
tions that are consistent considering all this multi-/} data. This 
approach is to be contrasted with previous modeling works that 
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Fig. 10. Range of allowed evolutions of the stellar mass den- 
sity computed with the range of cosmic SF histories of 
Fig. |5] The dashed lines corresponds to the best-fit to the IR 
counts. T he data points with error bars are measuremen ts from 
iPerez-Go nzalez et al. (2008) (em pty circles) and Eisner et al.l 
dlOOSl) (filled dots). The realistic (fialdrv & Glazebrookl ([2003h 
IMF is used both for the data points and to derive the stellar 
mass density from the SFR density. 



were based on predictions from a single preferred model: a range 
of models is given here. 

The inversion technique does not use any redshift informa- 
tion as input, although such an option can be (and has been) 
considered through a prior constructed on the low-redshift LF. 
Despite this arguably questionable lack of information about 
the redshift of the sources, the inversion technique recovers the 
known redshift distributions surprisingly well up to z = 2. The 
reason for this success probably lies in the very broad basis of 
wavelengths used in the inversion: 15;um, 160//m, and 850 /im 
(and intermediate wavelengths) do not probe the same popula- 
tions of galaxies at the same redshifts because of very different 
k-corrections. The uncertainties inherent in the library of SEDs 
that we used seem to be masked, at first order, by the extent of 
the data set that we considered. 

We find new constraints for the SFRD and its decomposition. 
Our method shows that the IR downsizing must be at work, even 
though only the IR counts are considered. Quantitatively, we are 
in good agreement with direct measurements of the SFRD at low 
and high redshifts. Again the clear advantage of our approach is 
its exhaustivity: the range of possible SF histories inferred from 
the inversion does not suffer from incompleteness, in contrast to 
surveys based on spectroscopic and photometric redshifts. This 
range matches recent measurements of the evolution of the stel- 
lar mass density very well, when a non-evolving IMF is used 
and stellar remnants are taken into account, in contrast to previ- 
ous works. 

Then by comparing the results on the SFRD from the "di- 
rect" (redshift-based) and "indirect" (counts inversion) methods, 
we showed that the population of HLIRGS tentatively seen at 
z > 2 is excluded, at least with SEDs similar to local ULIRGS. 
Either the photometric redshifts of these peculiar sources are 
systematically wrong, or their SED is very different from the 
local equivalents. Both reasons contribute probably to the ob- 
served discrepancy. 



Moreover, a strong contribution from obscured AGNs to the 
mid-IR counts, especially at lAjim as shown by toaddi etal] 
(2007a) cannot be excluded. For this reason, the objects that 
seem to be HLIRGS at z = 2-3 may also be actually ob- 
scured AGNs. Solid identifications of AGNs from far-IR diag- 
nostics will be needed for further modeling work. 

The set of galaxy SEDs that we used in the inversion (CEOl) 
was calibrated at z = for A < 100//m. Moreover, the SFR 
predictions made with t hese templates are compatible with radi o 
estimates up to z = 1 (lElbaz et al.ll2002t lAppleton etani2004l) . 
Although the multi-/l IR counts are reproduced reasonably well 
with the inverted luminosity functions, the 160 ;um counts cannot 
be satisfactorily reproduced by any SF history satisfying the rest 
of the counts. At this wavelength, the counts are systematically 
under-produced by models. As noted by LDP03, a population of 
objects colder than those included in the CEOl library is needed 
to produce these counts. We can therefore decide either on an 
evolution of the SEDs that were different - presumably colder 
for the same Lir - at high redshift, or on a poor calibration of 
the CEOl SEDs in this range (which is actually expected). This 
strongly confirms the tentative evidence fro m previous works 
on modeling (LDP03) or ob servation al Papov ich et al.l (l2007l) 
grounds. The recent works o f iPope et al, (.2008.) " on submillime- 
ter galaxies go in the same direction: these high-redshift IR 
galaxies are colder than previously thought. But only Herschel 
will enable the characterization of the dust temperature of these 
objects at these wavelengths. In fact, Herschel studies of local 
galaxies will be very helpful for high-redshift studies of galaxy 
evolution. 

Moreover, the inversion does not completely take the special 
case of AGNs into account. Implicitly, we assumed that their 
SEDs are close enough to the templates of star-forming galaxies 
represented in the CEOl library. Obviously, the AGNs have dif- 
ferent SEDs, with a flatter SED in the mid and far-IR. However, 
our results are not very sensitive to this issue. First, the contribu- 
tion of AGNs (as identified from X-rays or optical spectra) to the 
mid-IR counts is small, especially for faint fluxes. Therefore, a 
small correction to the bibliographic counts should be made, but 
remains difficult at all wavelengths. Moreover, as noted above, 
the global trends derived for the LF or the cosmic SFRD com- 
pare well with independent measurements at low redshift, which 
seems to indicate that the inversion is relatively insensitive to 
the precise shape of the SEDs. Of course, there are limits to this 
statement, but we believe that the small fraction of AGNs found 
in the counts, combined with this statement, should only weakly 
affect our findings. 

After checking that the predicted redshift distributions of the 
sources making the counts at 16, 24, and 70;um are compatible 
with the real redshift distributions, we then made predictions for 
the future surveys to come with Herschel, SPICA, SCUBA2, and 
Artemis (see Appendix ICl where we explain the bivariate distri- 
butions in z and Lir of the sources as a function of flux, and 
we give the fraction of the CIRB that will be resolved by future 
confusion limited surveys). 

Finally, we stress that, although the LFir and the SF history 
derived from the inversion actually depend (slightly) on the li- 
brary of SEDs chosen for the work, the use of a new, updated, 
SED library would be straightforward: unlike previous model- 
ing approaches, it would not imply a whole new work to refine 
parameters in order to obtain a good modeling of the counts. 
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Appendix A: The inverse problem 

A.1. The model 

As argued in the main text, (Sect. 2.1) the formal equation relat- 
ing the number of counts of galaxies N{Ai, S k) with the flux S u 
(within AS) at wavelength Aj (within dA) to the number of counts 
of galaxies, N{z, Lir), at redshift z (within dz) and IR luminosity 
LiR (within dLiR) is given by 
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N{A 



[Sk - F(Ai,z,LiR)]N(z,LiR)dzdLiR , (A.l) 



where 6d is the standard Dirac function, F the flux observed 
in a photometric band centered on wavelength Aj of a galaxy at 
redshift z with a Lir luminosity: 



F(Ai, z, Lir) = A/Dl(z)K(Ai, z, Ljr) . 



(A.2) 



Here, Di(z) is the luminosity distance for an object at redshift 
z with the standard cosmology used in this paper, A is the solid 
angle corresponding to one square degree, and K corresponds to 
the k-corrections: 



K(Ai,z,LiR) = l/R 



L^HA/(l+z)) 



1™.„ 1 + 



Ti{A)dA , 



(A.3) 



where Ti{A) is the transmission curve for the filter centered on Ai, 
R = Ti(A)dA, and L^'R(^) is the underlying hbrary of SEDs 

(CEOl) for which the SED of a galaxy depends only on its total 
luminosity Lir. 

As mentioned in the main text, from the point of view of the 
conditioning of the inverse problem, it is preferable to reformu- 
late Eq. (lA.U in terms of Z = logio(l + z), S = logio(5) and 
miR = logio Lir: 

N(Ai, Sk)^jj H(Sk, Ai, Z, m^mZ, miR)dZdmiR , (A.4) 

where the kernel of Eq. ( IA.4l i reads 

H(S,A,Z,myid = 10^-^'^+^+'"«5d[5 - F{A, 10^, lO'"")] 

with 

N{Z, myid = N{IQ^, 10'"™) , (A.5) 

N{Ai, Sk) = N(Ai, 10'^')10^-^'^' . (A.6) 

Here we have introduced the Euclidian-normalized number 
count, N, by multiplying the number count by the expected S 
power law. 

A.2. Discretization 

Let us project N(Z, '«ir) onto a complete basis of p x q functions 
{ei(Z)e/(miR))j=i,...,p , 

of finite (asymptotically zero) support, which are chosen here to 
be piecewise constant step functions: 



NiZ,miR) = ^ ^ ej(Z)eiimiR), 



(A.7) 



j=i 1=1 



The parameters to fit are the weights rtji. Calling 
- {'^ji)j=i,.-p,i=i,-q (the p X q parameters) and 
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Y = {N(Ai,Sii)]i=i_„r,k=i,..s (the r x s measurements), Eq. jAAj 
then becomes formally 

Y = M-X, (A.8) 
where M is a (r, s) x (p, q) matrix with entries given by 



M, 



ikjF 



I I ej(Z)eiimiR)HiSk,Ai,Z,miR)dZdmiA 
•J J ) ik jl 



A.3. Penalties 

Assuming that the noise in N can be approximated to be Normal, 
we can estimate the error between the measured counts and the 
non-parametric model by 



L(X) = x^C^) = (Y-M-X)^-W-(Y-M-X), 



(A.9) 



where the weight matrix W is the inverse of the covariance ma- 
trix of the data (which is diagonal for uncorrelated noise with 
diagonal elements equal to one over the data variance). Since we 
are interested here in a non-parametric inversion, the decomposi- 
tion in Eq. iA.lj typically involves many more parameters than 
constraints, such that each parameter controls the shape of the 
function, A^, only locally. As mentioned in the main text, some 
trade-off must therefore be found between the level of smooth- 
ness imposed on the solution in order to deal with the artefacts 
induced by the ill-conditioning, on the one hand, and the level 
of fluctuations consistent with the amount of information in the 
counts, on the other hand. Between two solutions yielding equiv- 
alent likelihood, the smoothest is chosen on the basis of the 
quadratic penalty: 



R(X) = X^-K-X, 



(A. 10) 



where K is a positive definite matrix, which is chosen so that 
R in Eq. ( lA.lOb should be non zero when X is strongly vary- 
ing as a function of i ts indices. In practice, we use a square 
Lapla cian pe nalization lOcvirk et al.l (D2 norm as defined by Eq. 
30 o f l2006bh . Indeed, a Tikhonov penalization does not explic- 
itly enforce smoothness of the solution, and a square gradient 
penalization favors flat solutions that are unphysical in our prob- 
lem. 

As mentioned in the main text, for a range of redshifts, a 
direct measurement, Xo, which can be used as a prior for A^, is 
available. We may therefore add as a supplementary constraint 
that 



P(X) = (X-Xo)^-W2-(X-Xo) 



(A.ll) 



should remain small, where the weight matrix, W2, is the inverse 
of the covariance matrix of the prior, Xq, and should be non zero 
over the appropriate redshift range. In short, the penalized non- 
parametric solution of Eq. dA.St accounting for both penalties is 
found by minimizing the so-called objective function 



Q{X)^L(X) + AR(X)+fiP{X), 



(A.12) 



where L(X), /?(X), and P(X) are the likelihood and regulariza- 
tion terms given by Eqs. ( |A.9t , dA.lOt , and ( lA.llt . respec- 
tively. The Lagrange multipliers A,jj > Q allow us to tune the 
level of smoothness of the solution (in practice, we set A = 0.02 
for the reasons given below) and the requirement that X should 
remain close to its prior for the range of redshifts for which data 
is available. The introduction of the Lagrange multipliers is for- 
mally justified by our wanting to minimize the objective function 



Q(X), subject to the constraint that L(X) and P(X) should faU in 

the range A^data ± V2 A^data and A^param ± V^^Vp^, respectively. 

The minimum of the objective function, QiX), given by 
Eq. ( IA.I2I 1 reads formally as 

X = (M-'-W-M+iK-HjuW2)"'-(M-'-W-Y-HjuW2-Xo) . 

This equation clearly shows that the solution tends towards Xq 
when yu ^ 00, while the smoothing Lagrange multiplier. A, 
damps counterparts of the compo nents of Y correspon ding to 
the higher singular vectors of M dOcvirk et al.l l2006bb . When 
dealing with noisy datasets, the non-parametric inversion tech- 
nique may produce negative coefficients for the reconstructed lu- 
minosity function. To avoid such effects, positivity is imp osed on 
those co efli cients «,/ in Eq. (|A.7| i, see for instance Ocvirk et al.l 
(l2006bl) or iPichon et all (|2002|) . In practice, the minimum of 
the objective function is found iteratively, using optimpack 
(Thiebaut 2005). The relative weight on the likelihood and the 
two penalties is chosen so that the three quantities have a compa- 
rable contribution to the total likelihood after convergence. This 
corresponds to a reasonably smooth variation in the LF both in 
redshift and Lir, and imposes a solution that is always within Icr 
of the observed low-redshift LF, Xq, when fj. is not set to zero in 
Eq. jAl2] l. 



Appendix B: Test of robustness 

To quantify the confidence level of the inversion technique, we 
test its robustness. Starting from an arbitrary LF, we produce 
IR counts in the bands and flux ranges corresponding to the ob- 
servations from this LF. Then, we add some random Gaussian 
noise to the simulated counts, using the real uncertainty on the 
observations as the cr of the error distribution for each flux bin. 
Finally, we apply the inversion technique described in Sect. 2.1 
to these noisy counts and obtain an output LF. 

The comparison of the input and output LFs is shown in 
Fig. IB. II The error on the absolute difference in logjo LFi„ - 
logjo LFout is represented in gray levels and contours. The differ- 
ence is generally less than 0.4 dex (factor 2.5) in the range where 
the LF can be constrained from the observed counts (range of the 
z-L plane encompassed by the dashed lines). A noticeable excep- 
tion is the very low-redshift range (z < 0.1), which corresponds 
to bright sources. For such large fluxes, the considerable noise 
in the observed counts produces large errors on the recovered 
LF. At high redshift, recall that the ultra-luminous population of 
galaxies appears as rare and very bright objects, in a flux range 
where the number counts are poorly known. 

Appendix C: IVIodel predictions for l-lerschel 

In Sect, in we have inverted the known IR counts to obtain con- 
straints on the evolving total IR LF. We have seen that the LF 
obtained through this inversion is realistic and matches most of 
the recent observations (counts, CIRB, Mid-IR LF at low red- 
shift). Then, in Sect. 15.21 we have shown how we can measure 
directly a part of this LF with a good confidence and that the 
LF resulting from the inversion is in good agreement with this 
solid measurement, validating the LF obtained by this empirical 
modeling approach. In this section, we use the median LF ob- 
tained from the inversions to predict some counts which should 
be observed with future observations in the FIR with Herschel 
or SCUBA2. 

At the time of publication, several new facilities are in prepa- 
ration to observe the Universe in the far-IR to sub-mm regimes. 
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Fig. B.l. Estimated robustness of the LF inversion used in this 
paper The relative difference between the input LF and the re- 
covered LF (when a realistic noise is added to the corresponding 
input counts) is larger for darker parts of the diagram. This differ- 
ence is relatively small (<0.4 dex) in the region of the z-L space 
effectively constrained by observations: the dotted and dashed 
lines correspond to the extreme fluxes considered at 24 yum and 
850//m, respectively, for this study. See main text for details. 
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Fig. C.l. Confusion limit for a 3.5 meter telescope. The 5cr- 
Ihour limits of SPICA-SAFARI (50//Jy, 33-210yum, dashed 
line), Herschel PACS (3mJy, 55-210yum, light blue line) and 
SPIRE (2 mJy, 200-670 /urn, blue line) are shown together with 
their wavelength ranges. The blue part of the curve is deter- 
mined by the source density criterion (i.e. the requirement to 
have less than 30% of the sources closer than O.SxFWHM), the 
red part is defined by the photometric criterion, i.e. sources must 
be brighter than 5 times the rms due to very faint sources below 
the detection limit. 



The definition of the confusion limit is not trivial, in particu- 
lar because it depends on the level of clustering of galaxies; the 
optimum way to define it would be to perform simulations to 
compute the photometric error as a function of flux density, and 



The differential counts (normalized to Euclidean) at wavelengths 
ranging from 16 to 850 /zm, which we derived from the inver- 
sion technique, are presented in Fig. IC.3I The separation of the 
contribution of local, intermediate, and distant galaxies in differ- 
ent colors illustrates the expected trend that larger wavelengths 
are sensitive to higher redshifts, hence the relative complemen- 
tarity of all IR wavelengths. There will be a bias towards more 
luminous and distant objects with increasing wavelength, illus- 
trated here for the Herschel passbands (see Fig. IC.4l i. but this 
may be used to pre-select the most distant candidates expected 
to be detected only at the largest wavelengths. In the following, 
we discuss the predictions of the inversion technique for those 
instruments, as well as their respective confusion limits, which 
is the main limitation of far-IR extragalactic surveys. 

The ESA satellite Herschel is scheduled to be launched 
within the next year, while the next-generation IR astronomi- 
cal satellite of the Japanese space agency, SPICA, is scheduled 
for 2010, with a contribution by ESA under discussion, includ- 
ing a mid-IR imager named SAFARI. Both telescopes share the 
same diameter of 3.5 meters, but the lower telescope tempera- 
ture of SPICA, combined with projected competitive sensitiv- 
ities, will make it possible to reach confusion around 70;um 
(where Herschel is limited by integration time). The 5cr-lhour 
limits of the instruments SAFARI onboard SPICA (50/iJy, 33- 
210yum, dashed line), PACS (3mJy, 55-210yL(m, light blue line) 
and SPIRE (2 mJy, 200-670/im, blue line) onboard Herschel are 
compared in Fig. lC.ll to the confusion limits that we derive from 
the best-fit model of the inversion, at all wavelengths between 30 
and 850 ;um, assuming the the confusion limit definition given 
below. 
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Fig. C.l. Detection limits for confusion limited surveys from 70 
to 850//m. The curves show the minimum IR luminosity (8- 
1000 yum), or equivalently SFR (=Lir x 1.72 10""'), that can be 
detected for a star-forming galaxy assuming that it has an SED 
similar to the Chary & Elbaz (2001) ones. The 70, 100, 160, 250, 
350 and 500 fim limits correspond to a 3.5 m telescope diameter, 
such as Herschel or SPICA, while the 400 /im is for a 12 m class 
telescope such as APEX (e.g. ARTEMIS, we show the average 
between the two bands at 350 or 450 fim to avoid confusion with 
Herschel) and the 850jum is for a 15 m telescope as the JCMT 
(SCUBA). 
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then decide that the confusion limit is e.g. the depth above which 
68 % of the detected sources are measured with a photometric 
accuracy better than 20 %. In the following, we only consider 
a simpler approach that involve s computing the tw o sources of 
confusion that were discussed in iDoleet al.l(l2003h : 

- the photometric confusion noise: the noise produced by 
sources fainter than the detection threshold. The photometric 
criterion corresponds to the requirement that sources are 
detected with an S/N(photometric)>5. 

- the fraction of blended sources: a requirement for the quality 
of the catalog of sources will be that less than N% of the 
sources are closer than O.SxFWHM, i.e. close enough to not be 
separated. 

We tested various levels for and found that N was 
equivalent to the above requirement that 68 % of the detected 
sources are measured with a photometric accuracy better than 
20 % using realistic simulations in the far IR for Herschel. We 
therefore use the value =30 %. The confusion limit is then de- 
fined as the flux density above which both criteria are respected. 
As a result, it is found that the main limitation is the fraction 
of blended sources at A-5Q-\Q5 jiva (blue part of the curve in 
Fig. IC. Il l and the photometric confusion noise below and above 
this range, i.e. at A= 33-50 and /l=105-210jum (red parts of curve 
in Fig. IC.lb . As a result of their smaller beam, shorter IR wave- 
length are more efficient at detecting faint star-forming galaxies 
than longer ones (see Fig. |C.2| i. This is at the expense of observ- 
ing farther away from the peak of the far IR emission, which 
implies larger uncertainties on the derivation of the total IR lu- 
minosity due to the uncertain dust temperature. 

We note that the confusion limit for a 3.5m-class telescope, 
such as Herschel, is ten times more than the depth it can reach in 
one hour (5cr). With a source density of 12.8 sources per square 
degree at the 500 /vm confusion limit (25 mJy), or equivalently 
one source in a field of 17 arcmin on a side, this shows that 
the best strategy at this wavelength is to go for very large and 
moderately shallow surveys, in order to identify the rare and very 
luminous distant galaxies. 

For comparison, we also illustrated the ground-based ca- 
pacity of ARTEMIS built by CEA-Saclay which will operate 
at the ESO 12 m-telescope facility APEX (Atacama Pathfinder 
Experiment) at 200, 350 and 450yL/m and SCUBA-2 that will 
operate at the 15 m telescope JCMT at 450 and 850yum. To 
avoid confusion between all instruments, we only show the av- 
erage wavelength 400 jum for a 12 m-class telescope and 850 yum 
for a 15 m-class telescope (Fig. IC.lt . Although the confusion 
limit in the 850 ;um passband is ten times below that of Herschel 
at the longest wavelengths, this band is not competitive with 
the ~400 jum one, which should be priorities for ARTEMIS and 
SCUBA-2 for the study of distant galaxies, or with the 70 and 
100 //m ones for a 3.5 m space experiment such as SPICA and 
Herschel, for redshifts below z ~5. We also note that only in 
these two passbands will the cosmic IR background be resolved 
with these future experiments (see Table ICTT l. which suggests 
that a larger telescope size should be considered for a future ex- 
periment to observe the far IR Universe above lOOjum and be- 
low the wavelength domain of ALMA. We did not mention here 
ALMA since it will not be affected by these confusion issue: 
due to its very good spatial resolution, it will be limited to either 
small ultradeep survey, hence missing rare objects or follow-ups 
of fields observed with single dish instruments, e.g. ARTEMIS. 
Finally, the JWST that will operate in the mid IR will be a very 
powerful instrument for probing the faintest star-forming galax- 
ies in the distant Universe, but predictions are difficult to produce 



at the present stage since it has already been found that extrapo- 
lations from the mid to far IR become less robust already at z ~1 
(e.g. ,Daddi et aL.2007k.Papovich et al...2007;,Pope et aL.2008l) . 
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Table C.l. Fraction of the CIRB resolved by confusion-limited Herschel surveys 
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Fig. C.3. Counts predicted from the inversion in the far-infrared and sub-mm (solid line). The counts are decomposed in redshift 
bins (blue=z < 0.5; green=0.5 < z < 1.5; orange=1.5 < z < 2.5; red=z > 2.5). The oblique dashed line coiTesponds to the limit in 
statistics due to the smallness of a field like GOODS North-nSouth or 0.07 square degrees: less than 2 galaxies per flux bin of width 
(5 logF=0. 1 dex are expected below this limit. 
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Fig. C.4. Differential counts predicted from the non-parametric inversion for future Herschel observations: PACS 100//m (soHd 
black), SPIRE 250 fim (dotted blue), SPIRE 350 jum (dashed green), SPIRE 500 yum (dot-dashed red) decomposed simultaneously 
in redshift and Lir. 



